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Abstract 

We propose a novel method for constructing Hilbert transform (HT) 
pairs of wavelet bases based on a fundamental approximation-theoretic 
characterization of scaling functions — the B-spline factorization theorem. 
In particular, starting from well-localized scaling functions, we construct 
HT pairs of biorthogonal wavelet bases of L 2 (R) by relating the corre- 
sponding wavelet filters via a discrete form of the continuous HT filter. 
As a concrete application of this methodology, we identify HT pairs of 
spline wavelets of a specific flavor, which are then combined to realize a 
family of complex wavelets that resemble the optimally-localized Gabor 
function for sufficiently large orders. 

Analytic wavelets, derived from the complexification of HT wavelet 
pairs, exhibit a one-sided spectrum. Based on the tensor-product of such 
analytic wavelets, and, in effect, by appropriately combining four separa- 
ble biorthogonal wavelet bases of L 2 (E 2 ), we then discuss a methodology 
for constructing 2D directional-selective complex wavelets. In particu- 
lar, analogous to the HT correspondence between the components of the 
ID counterpart, we relate the real and imaginary components of these 
complex wavelets using a multi-dimensional extension of the HT — the di- 
rectional HT. Next, we construct a family of complex spline wavelets that 
resemble the directional Gabor functions proposed by Daugman. Finally, 
we present an efficient FFT-based filterbank algorithm for implementing 
the associated complex wavelet transform. 

1 INTRODUCTION 

The dual-tree complex wavelet transform (DT-CWT) is a recent enhancement of 
the conventional discrete wavelet transform (DWT) that has gained increasing 
popularity as a signal processing tool. The transform was originally introduced 
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by Kingsbury [17, 19] to circumvent the shift-variance of the decimated DWT, 
and involved two DWT channels in parallel with the corresponding wavelets 
forming a quadrature pair. In particular, Kingsbury realized the quadrature 
relation by interpolating the lowpass filters of one DWT "mid-way" between the 
lowpass filters of the other DWT. Moreover, based on appropriate combinations 
of separable wavelets, he extended the dual-tree construction to two-dimensions, 
where the corresponding transform, besides improving on the shift-invariance of 
the 2D DWT, exhibits better direction selectivity as well. There is now good 
evidence that the transform tends to perform better than its real counterpart 
in a variety of applications such as such as deconvolution [10], denoising [29], 
and texture analysis [16]. 

The crucial observation that the dual-tree wavelets involved in Kingsbury's 
construction form an approximate HT pair was made by Selesnick [28, 26]. He 
also demonstrated that a particular phase relation between the lowpass (refine- 
ment) filters of the two channels resulted in the desired HT correspondence. 
This link consequently transposed the problem of designing different flavors of 
dual-tree wavelets to that of identifying new HT pairs of wavelets. Indeed, fol- 
lowing this remarkable connection, several new paradigms and extensions have 
been proposed: design of HT pairs of biorthogonal wavelet bases [40], alterna- 
tive frameworks for complex non-redundant transforms [12], and the M-band 
extension [8], to name a few. 

1.1 Motivation 

The deployment of complex signal representations for the determination of in- 
stantaneous amplitude and frequency is classical [13, 38]. Gabor and Villc 
[13, 39] proposed to unambiguously define them using the concept of the an- 
alytic signal — a unique complex-valued signal representation specified using 
the HT. Specifically, the analytic signal s a (x) = s(x) + jHs(x) correspond- 
ing to a real-valued signal s(x) (H denotes the HT operator), was used to 
stipulate the instantaneous amplitude and phase via the polar representation 
s a (x) = A{x)e^^ . In particular, this representation allows one to retrieve 
the time-varying amplitude and frequency of an AM-FM signal of the form 
s(x) = A(x)cos(2n J x v(f)dT + £o) via the estimates ^4 est (x) = |s a (ar)| and 
f est (x) = (27r) _1 d(p{x) I dx , assuming A(x) to be slowly-varying compared to 
v(x). The analytic signal has become an important complex- valued representa- 
tion in signal processing, especially in applications such as phase and frequency 
modulation, speech recognition and processing of seismic data. These concepts 
have also been transposed to the multi-dimensional setting: the local frequency 
has been used as a measure of local signal scale; structures such as lines and 
edges have been distinguished using the local phase; and the local amplitude 
and phase have been used for edge detection and for texture and fingerprint 
analysis [25]. 

The advantage of viewing the dual-tree wavelets as a HT pair is that we can 
make a direct connection with the formalism of analytic signals. Indeed, if we 
transpose the above concept to the wavelet domain and consider the input signal 
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to be locally of the AM-FM form, we obtain a response where the local energy 
of the signal is encoded in the magnitude of the wavelet coefficients, while the 
relative displacement is captured by the phase. In fact, this turns out to be 
the fundamental reason for the superiority of the DT-CWT over conventional 
real-valued transforms whose response is necessarily oscillating. 

1.2 Our Contribution 

In this contribution, we invoke the B-spline factorization theorem [37] — a fun- 
damental spectral factorization result — along with certain fractional B-spline 
calculus [36] , to construct HT pairs of biorthogonal wavelets from well-localized 
scaling functions. In particular, we do so by relating the corresponding wavelets 
filters via a discrete version of the continuous HT filter. 

Next, we identify a family of analytic spline wavelets, of increasing vanishing 
moments and regularity, that asymptotically converge to Gabor-like functions 
[13]. As far as the implementation is concerned, unlike Kingsbury's scheme that 
uses different filters for different stages (often with filter-swapping between the 
dual-trees), our implementation uses the same set of filters at all stages of the 
filtcrbank decomposition. Notably, we use an appropriate pair of projection 
filters for coherent signal analysis which, in turn, allows us to identify a discrete 
counterpart of the analytic wavelet — the so-called analytic wavelet filter that 
exhibits a one-sided spectrum. 

The construction is then extended to two-dimensions through appropriate 
tensor-products of the one-dimensional analytic wavelets. In particular, we con- 
struct a family of directional complex wavelets that resemble the directional 
Gabor functions proposed by Daugman [9] for sufficiently large orders. More- 
over, we also relate the real and imaginary components of the complex wavelets 
using the directional HT — a multidimensional extension of the HT — that pro- 
vides further insight into the directional-selectivity of the dual-tree wavelets. 

1.3 Organization of the Paper 

We begin by recalling certain fundamental definitions and properties pertain- 
ing to the HT and the fractional B-splines in §2. We characterize the action 
of the HT operator on B-splines in §3, which, along with the B-spline factor- 
ization theorem, is used to propose a formalism for constructing HT pairs of 
biorthogonal wavelet bases in §4. The implementation aspects are discussed 
in §5. As a concrete application, we construct the Gabor-like wavelets in §6. 
In §7, directional complex wavelets are constructed by appropriately combining 
the wavelets corresponding to certain separable multiresolution analyses; the 
highlight of this section is the construction of 2D Gabor-like spline wavelets. 
The implementation aspects of the corresponding 2D Gabor-like transform are 
provided in §8, before concluding with §9. 
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2 PRELIMINARIES 



We begin by introducing specific operators and functions that play a major 
role in the sequel followed by a discussion of their relevant properties. In what 
follows, we use /(<*;) = J Rd f(x)e~^ x w dx to denote the Fourier transform of 
a function f(x) on R d (d ^ 1), with x T u being the usual inner-product on R d . 
We also frequently use the notations /(• — s) and /(A-), corresponding to some s 
in R d and A > 0, to denote the function obtained by translating (resp. dilating) 
f(x) by s (resp. A). We denote the Kronecker-delta sequence by S[n]: its value 
is 1 at n = 0, and is zero at all other integers. 

2.1 Hilbert Transform and Wavelets 

The Hilbert transform, that generalizes the notion of the quadrature transfor- 
mation cos(luqx) sin(wox) beyond pure sinusoids [5], forms the cornerstone 
of this paper. From a signal-processing perspective, the HT can be interpreted 
as a filtering operation in which the amplitude of the frequency components is 
left unchanged, while their phase is altered by ±7r/2 depending on the sign of 
the frequency. 

Mathematically, the HT of a sufficiently well-behaved function is defined 
using a singular integral transform [2, 30]. However, in the context of finite- 
energy signals, it admits a particularly straightforward formulation based on 
the Fourier transform on L 2 (M). In particular, the Hilbert transform on L 2 (R) 
is characterized by the equivalence 

Hf(x) ^ -jsignM/M (1) 
where the multiplier sign(w) is defined as for non-zero oj, and as zero at 

LO = 0. 

Based on the above definition 1 , and the properties of the Fourier transform 
on L 2 (K), the following properties of the HT can be readily derived: 

• Linearity and Translation-invariance: It is a linear and translation-invariant 
operator; that is, it acts as a convolution operator. 

• Dilation-Invariance: It commutes with dilations: H{.f(X-)}(x) = (Hf)(\x)for 
all A > 0. 

• Anti- Symmetry: It anti-commutes with the flip operation f T (x) = f(—x), 
so that (Hf T )(x) = —(Hf) T (x); thus the HT of a symmetric function is 
necessarily anti-symmetric. 

• Unitary (Isometric) Nature: It acts as a unitary operator on L 2 (M), so 
that (Hf,Hg) — (f,g) for all / and g, where (•, •) denotes the usual inner- 
product on L 2 (M). Equivalently, this means that the inverse HT operator 
is given by its adjoint: H^ 1 = H*. 

1 The definition can also be extended to tempered distributions such as the Dirac delta and 
the sinusoid [2, §2.5]. 
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It is well-known that HT of a wavelet is also a wavelet. The implication 
of the simultaneous invariance to dilations and translations is that the HT of 
a dilated-translated wavelet is a wavelet, dilated and translated by the same 
amount: H{ip(\x — s)} = (Ti.ip)(Xx — s). Moreover, an immediate consequence 
of the unitary property is that the HT operator maps a basis into a basis: if {ip n } 
form a (Riesz) wavelet basis of L 2 (R), then so does {Hip n }. It even preserves 
biorthogonal wavelet bases of L 2 (M): if {ip n } and {ip m } form a biorthogonal 
wavelet basis of L 2 (M), satisfying the duality criteria (ip n , ip m ) = S[m — n], then 
using the same unitary property, we have 

{Hip n ,H^ m ) = (ipn,i>m) = S[m- n] (2) 

so that {Hipn} and {Htp m } form a biorthogonal wavelet basis of L 2 (M) as well. 
It is exactly the above invariance properties that make the construction of HT 
pair of wavelet bases of L 2 (R) feasible. 

Unfortunately, the HT exhibits certain inherent pathologies in the context 
of multircsolution analyses and wavelets. The impulse response of the HT, 
TL8(x) = 1 / ttx (in the sense of distributions) , clearly indicates the non-local 
nature of the operator. This has two serious implications: (i) the HT of a 
compactly-supported scaling function/ wavelet is no longer of finite support; (ii) 
the HT-transformcd function has a 0(l/|a:|)-decay in general, and hence is not 
intcgrablc; and (Hi) the (anti-symmetric) HT suppresses the dc-component of 
symmetric scaling functions that is essential for fulfilling the partition-of-unity 
criterion. Therefore, the HT of a scaling function is not a valid scaling function, 
and cannot be used to specify a multiresolution analysis in the sense of Mallat 
and Meyer [22, 23]. 

Next, we recall the notion of an analytic signal that generalizes the phasor 
transformation transformation cos(u)ox) i— > cxp(ja;o2;) to finite-energy signals 
using the HT as the quadrature transformation. In general, the analytic signal 
fa{x) associated with a real- valued signal f(x) is defined as the complex- valued 
signal 

f a (x) = f(x) + jHf(x). (3) 

In particular, f a (x) = exp(jiv x) when f(x) — cos(uoqx). Importantly, note that 
the Fourier transform of the analytic signal evaluates to f a (oj) = (1 + sign(w)) f(u>), 
so that f a (uj) vanishes for all negative frequencies. It is exactly this one-sided 
spectrum that makes the analytic signal particularly interesting in signal pro- 
cessing [13]; we exploit this property for constructing directional wavelets in 
§7- 

2.2 Fractional B-spline Multiresolution 

The family of fractional B-splines [4] — fractional extensions of the polynomial 
B- splines — will play a key role in the sequel. In particular, we recall that the 
fractional B-spline (3"(x), corresponding to a degree a € Kg an< ^ a snn?t T *= ^> 
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is specified by its Fourier transform 

The parameters a and r control the width and the average group delay of the 
scaling function respectively. In particular, when r = (a + l)/2, the fractional 
B-spline corresponds to the causal B-spline /3+(x) defined in [36]. The 

fractional B-splines, in general, do not have a compact support (except for 
integer degrees); however, their 0(l/|:r| Q+2 ) decay ensures their inclusion in 
L 1 (M) n L 2 (M). Another relevant property that will be invoked frequently is 
that the shift r influences only the phase of the Fourier transform; that is, 
1/3" (w) | is independent of r. 

The fundamental role played by fractional B-splines in this paper is, however, 
based on the fact they satisfy certain admissibility criteria [36, 4] needed to 
generate a valid multiresolution of L 2 (M): 

(CI) The approximation space Vo = span £2 {/3"(- — fc)}fcez admits a stable Ricsz 
basis. 

(C2) There exists an intcgrable sequence h"[k] (refinement filter) such that the 
two-scale relation 

\pr(D=Y. h r[m{*-k) (5) 

fcez 

holds. In particular, the transfer function of the refinement filter is specified by 

H?(en = ^ (1 + e^)^- T (1 +c-^)^ +T . (6) 

(C3) Partition of unity: The integer-translates of (i"{x) can reproduce the unity 
function. 



We briefly discuss the significance of these admissibility conditions. The 
criterion (CI) ensures a stable and unique representation of functions in Vo 
using coefficients from I 2 {It); cquivalcntly, this also signifies that the transfer 
function of the autocorrelation (Gram) filter, A a (e : ' UJ ) = X^feez \Pt{ IjJ + 27rfc)| 2 , 
is uniformly bounded from above, and away from zero [21]. On the other hand, 
(C2) implies the inclusion of [3^(x/2) in Vo, which, in turn, allows one to define 
a hierarchical embedding of approximation spaces {Vj}j<£z that is key to the 
multiresolution structure of the associated wavelet transform. Finally, the tech- 
nical condition (C3) ensures that the multiresolution {V,} is dense in L 2 (M): 
arbitrarily close approximations of functions in L 2 (M) can be achieved using 
elements from {V,}. 



3 HILBERT TRANSFORM AND B-SPLINES 

It turns out that the action of the HT on B-splincs can be effectively character- 
ized in terms of certain fractional finite-difference (FD) operators. In particular, 
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corresponding to an order a G Rq" and shift t € K, we consider the operator 
A" defined on L 2 (K) by 

A«/(x) ^ D?(ei")f(w), (7) 



where D?(e joj ) = (l - c"^) 2 (l 

One recovers the conventional n-th order FD operator by setting a = n 
and t = n/2. Since the operator has a periodic frequency response, one can 
associate with it a digital filter d"[k] through the correspondence A"f(x) = 
d"[k]f(x — k). The FD operator that is especially relevant for our purpose is 
the zcroth-order operator A°_ (henceforth, we simply denote it by A). The 

corresponding frequency response D(e JW ) = Z?°_ 1/ . 2 (e : ' w ) reduces 2 to 

D(e> u ) = -jsign(w)e"^ for lu e (-tt, tt), (8) 

signifying that D{e^) is in L 2 ((— 7r, 7r)); the corresponding filter coefficients 
d[fc] = d° 1/2 W m ^ 2 (^) are then specified 3 by 

d[k) = ^- f D{e^)e^duj 

J — 7T 

Thus, similar to the HT operator, A is also unitary, and the corresponding filter 
d[k] can be interpreted as a discrete form of the continuous HT filter 1/ttx. In 
particular, we can relate the action of the HT on the B-splines solely in terms 
of this filter. Indeed, it can easily be seen that the Fourier transform of the 
B-spline can be factorized as 

#?M = U^) 1/2 (-^)- 1/2 D(el")^ +1/2 (w), (10) 

which, along with the identity (jui)i (—joj)~i = jsign(w), results in the equiv- 
alence 

Hff{x) -jsign(w)^(a;) 

= -3 sign(w) • jsign(w)L>(e^)/^ +1/2 (w) 

= D(end? +1/2 n 

^A/3« +1/2 (x), (11) 
that establishes the desired result: 



2 We specify the fractional power of a complex number z by z~> = |z| 7 e J7 ar s( z ) corre- 
sponding to the principal argument |arg(z)| < tt. On this principal branch, the identity 
(21Z2) 7 = zjz^ holds only if arg(zi) + arg(zi) 6 (— it, it) [31, Chapter 3]. 

3 The inverse Fourier transform over the principal period (—tt, tt) is invoked. 
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Proposition 3.1 The HT of a fractional B-spline can be expressed as 



W*) = E ^kTW)^ +1/2(x k) ' (12) 

In particular, the digital filter d[k] acts as a unitary convolution operator on 
L 2 (]R) when applied to functions, and as a discrete filter on ^ 2 (Z) when applied 
to sequences. The theoretical difficulty with the HT stems from the fact that 
its frequency response has a singularity at ui = 0, which results in a poor decay 
of the transformed output. The remarkable feature of (12) is that we have 
been able to express the slowly decaying HT as a linear combination of the 
better-behaved B-splines. Specifically, the sequence d[k] decays only as 0(l/\k\), 
whereas /3" +1 / 2 ( x ) decays as 0(l/\x\ a+2 ). 

Thus, by expressing the HT using shifted B-splines as in (12), we have, in 
effect, moved the singularity onto the digital filter. In the sequel, we shall apply 
this filter to the wavelets where its effect is much more innocuous since V>(w) = 
around the origin. 

Half-Delay Filters: As remarked earlier, the shift parameter r only affects 
the phase of the Fourier transform of the fractional B-spline and the correspond- 
ing refinement filter [4]. In particular, based on the factorization 

= (! +e^)^ Mr+ ^ ) (l + e -^)^ +(r+i) 
= (1 + e J ' w )- 1 / 2 (l + e-^y^H^i^) 
= e- j %H?(e j "), for cj e (-tt.tt), 

we arrive at the following result: 

Proposition 3.2 The spline refinement filters h"[k] and h" +1 ^ 2 [k] are "half- 
sample" shifted versions of one another in the sense that 

H? +1/2 (en=e- j ^H?(en (13) 

for all oj in (—it, it). 

Indeed, if we consider the bandlimited function h"(x) — h"[k]s'mc(x — k) 
that satisfies the constraint h"(x)\ x= k = h"[k], then we have, as a consequence 
of (13), the relation h" +1 , 2 [k] = h"(k— 1/2): each filter provides the bandlimited 
interpolation of the other mid-way between its samples. 

Finally, we make a note of the fact that the above refinement filters can also 
be related through a conjugate-mirrored version of the FD filter: 

H? +1/2 (en = D(-e~nH?(en- (14) 
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4 HT PAIR OF WAVELET BASES 



Before stating the main results, we recall the approximation-theoretic notion of 
approximation order, and a fundamental spectral factorization result involving 
B-splines. 

Approximation Order: Scaling functions play a fundamental role in wavelet 
theory. The technical criteria for a valid scaling function was discussed earlier 
in the context of B-splincs (cf. §2.2). Next we recall the fundamental notion 
of order for a scaling function that characterizes its approximation power [37]. 
A scaling function <p(x) is said to have an approximation order 7 if and only 
if there exists a positive constant C such that for all elements of the Sobolev 
space W2QR), of order 7, we have the estimate 

\\f-P a f\\^Ca^Wf\\. (15) 

Here P a denotes the projection operator from W^QR) onto the approximation 
subspace span^2{cp(-/a — k)}k&, and d 1 denotes the (distributional) derivative 
of order 7. In other words, the approximation order provides a characterization 
of the rate of decay of the approximation error for sufficiently regular functions 
as a function of the scale. 

It turns out that, akin to their polynomial counterparts, the order of frac- 
tional B-splines is entirely controlled by their degree [36, 4]; in particular, we 
have 7 = a + l. Equivalently, this signifies that any polynomial of degree ^ \a~\ 
can be reproduced by the set {/3"(- — k)}, which is crucial for capturing the 
lowpass information in images is concerned. 

Characterizazion of Scaling Functions : A fundamental result in wavelet the- 
ory is that it is always possible to express a valid scaling function as a convolu- 
tion between an fractional B-spline and a distribution [37] . The original result 
in [37] involves causal B-splines; however, the result can readily be extended 
to the more general fractional B-splines since the shift parameter r docs not 
influence the order of the scaling function. Indeed, note the theorem in [37] 
asserts that H(e : ' UJ ) is the refinement filter of a valid scaling function (cf. §2.2) 
of order a + 1 if and only if it can be factorized as 

/3" spline part 
^ s\ ^distributional part 

ff(0=(^— ) Q(en , (16) 

where Q{e ju) ) is stable: |Q(e iw )| < C < +00 for all w. Rewriting (16) in 
terms of a (a, r) B-spline refinement filter, we then have the following equivalent 
representation: 

/3° spline part 

/ a* ^distributional part 

#(en=( i± |— ) P(en (17) 
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with P(e ju ) = e"^(^" T )g(e^)for u E (-7r,7r). Note that P(e^) is stable, 
with |-P(e iw )| < C < +00 for all uj. That is, H(e lu> ) is the refinement filter of a 
valid scaling function of order a + 1 if and only if it admits a stable factorization 
as in (17). We then arrives at the following extension: 

Theorem 4.1 fB-spline Factorization ) A valid scaling function cp(x) is of or- 
der a + 1 if and only if its Fourier transform can be factorized as 

?H = ^(u)$oM (18) 

for some r6l, where cpo(w) * s a function ofui that is bounded on every compact 
interval, and equals unity at the origin. 

In the signal domain, this corresponds to a well-defined convolution cp(x) = 
(/?" * cpo)(x) between a B-spline and the tempered distribution cp . The crux of 
the above result is that it is the constituent B-spline that is solely responsible for 
the approximation property, and other desirable features of the scaling function 
[37]. 

4.1 Construction of HT Pairs of Wavelets 

In what follows, we use the notation /^(x), corresponding to a function /(x), 
and integers j and k, to denote the (normalized) dilated-translated function 
2J/ 2 /(2 J x — k). The HT of a wavelet is also a wavelet in a well-defined sense. In 
particular, if tp(x) is a wavelet whose dilations-translations {ipj t k} form a Riesz 
basis of L 2 (M), then H.ip(x) is also a valid wavelet with {l~tipj,k} constituting a 
Riesz basis of L 2 (]R). As remarked earlier, this follows from the fundamental 
invariance properties of the HT. 

We now establish a formalism for constructing the HT of a given wavelet 
ip{x). In particular, if cp(x) be the associated scaling function, say of order a+1, 
and g[k] be the generating wavelet filter, then we have the relation 

1>(x/2) = Y,9[kMx-k). (19) 

Following Theorem 4.1, let us factorize cp(x) as 

cp(x) - (/?« * cp )(x) (20) 

corresponding to some real r. Then, consider the scaling function cp'(x), of the 
same order, specified by cp'(x) = (/9° +1 / 2 * ( Po)(x). Let ip'(x) be any arbitrary 
wavelet, corresponding to the multiresolution analysis associated with cp'(x), 
that is specified by ip'(x/2) = ^2 g'[k](p'(x — k). We then have the following 
necessary and sufficient condition for the desired HT correspondence in terms 
of the discrete HT filter d[k] (see §11.1 for a proof): 

Theorem 4.2 (WF Pair of Wavelets,) The wavelets ip(x) and tp'(x) have the 
correspondence tp'(x) — Hip(x) if and only if g'[k] = (d* g)[k]. 
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Moreover, the construction has the following characteristics: 

• Both cp(x) and cp'(x) have the same Riesz bounds and the same decay, 

• The refinement filters H(e' u ) and H'fe^) corresponding to (p(x) and 
cp'(x) respectively, are related as 

for all lo in (— tt,tt). 

The equality of the Riesz bounds follows from the observation that the autocor- 
relation filters of cp(x) and cp'(x) are identical. Indeed, we have 

a[k] = <cp, cp(- - k)) = (f3 2 a+1 * cp * cpo )(*), 
o'[fc] - (cp', <p'(- - k)) = (/3 2Q+1 * cpo * cpj)(fc). 

The assertion regarding the decay is based on the observation that both (3"(x) 
and /3" +1 / 2 ( x ) navc the same decay Finally, using (13) and (17), we can relate 
the transfer functions on (— ir, n) as follows 

H'(en = H? +1/2 (ei»)P(en 

where P(e 3 ") denotes the transfer function of the filter associated with the 
distribution cp . 

Remark: Note that although Hip(x) is unique, the scaling function cp'(x) 
and the corresponding filter g'[k] generating Hip(x) are by no means unique. 
For instance, the particular choice cp'(x) = H(p(x) and g' = g is sufficient to 
ensure that ip'(x) — Hip(x). Moreover, if <p'(x) and g'[k] generate the wavelet 
ip'(x/2) = ^g'[k\(p'{x — k) such that tp'(x) = Hip(x), then so do (p' eq (x) = 
J2r[k](p'(x — k) and g' cq [k] — (g' * ri nv )[fc]. Here the filter r[k] is such that 
< | J2k r [k]e~ : > u ' k \ < +oo for all lo so that the convolutional inverse n nv [fc] is 
well-defined. 

The condition g' [k] = (d* g) [k] is both necessary and sufficient only for our 
preferred choice of the scaling function cp'(x) = (/3" +1 ^ 2 *(p )(x). This particular 
choice of the scaling function against the more direct choice Tt(p(x) is justified 
on the following grounds: 

• The function cp'(x) is well-localized with better decay properties than 
Wcp(x); the latter is not even intcgrablc in general (e.g., the Harr scaling 
function) , 

• The scaling function (p'(x) satisfies the partition-of-unity requirement, 
whereas H<p(x) is not a valid scaling function since H(p(0) is not nec- 
essarily unity. For example, if cp(x) is symmetric and Hcp(x) is intcgrablc, 
then we have H(p(Q) = f(Hq))(x)dx = following the fact that Hcp(x) is 
anti-symmetric. 
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4.2 HT Pairs of Biorthogonal Wavelets 

A biorthogonal wavelet basis of L 2 (R) , corresponding to the dual-primal scaling 
function pair (cp,cp) of order (N + 1,N +1), involves the nested multiresolution 

{0} C • • • C V_i C V C Vi C • • • C L 2 (K), 

and its dual 

{0} C • • • C V-i C V C Vi C • • • C L 2 (M), 

where the approximation subspace Vj (resp. Vj) is generated by the translations 
of (pjfl(x) (resp. q)jfi(x)) [21]. Let (ip,ip) be the wavelets associated with 
these multiresolutions, which, along with their dilated-translated copies, encode 
the residual signal — the difference of the signal approximations in successive 
subspaccs. In particular, the wavelet tpjfi(x) (resp. Vj,o(x)) and its translates 
span the complementary space Wj = Vj Q Vj-i (resp. Wj = Vj Vj-i)- The 
crucial aspect of the construction is that the dilated-translated ensemble ipj : k(x) 
and ipj',k'{x) form a dual basis of L 2 (R), i.e., they satisfy the biorthogonality 
criteria {ipj,k>i ) j',k') = &\j — j\ k — k']. The expansion of a finite-energy signal 
f(x) in terms of this biorthogonal basis is then given by 

f{x)= (f>hM,k(*)- (21) 

(j,fe)ez2 

In other words, the wavelets k{ x )} an d {ipj,k( x )}, interpreted as the analysis 
and synthesis wavelets respectively, together constitute a biorthognal wavelet 
basis of L 2 (M). 

In particular, let <p(x) and cp(x) be the scaling functions, of order N + 1 
and N+l respectively, associated with a given biorthogonal wavelet basis, with 
associated wavelets 

i>{x/2) = J29[k](p{x-k), ^{x/2) = Y,9[kMx-k). 

Now, let cp(x) = (/3f * <f>o)( x ) an d <p{x) = (P^ * <Po){x) be the respective fac- 
torizations of <p(x) and cp(x). Consider the scaling functions cp'(x) = (/3f +1 / 2 * 
(f>o)(x) and (p'(x) = {fl^ +l / 2 * ( Po)( x ), with associated wavelets specified by 

V/(x/2) = E 9'[kW(x - k), V'(x/2) = 9'[kW(x - k). 

Then the following result comes as a direct consequence of Theorem (4.2). 

Corollary 4.3 (W£ Pair of Biorthogonal Waveletsj The following are equiva- 
lent: 

• The primal and dual wavelets form HT pairs, V>'(x) = Htp(x) and tp'(x) = 
TLip{x), and {ip'j k (x)} and {ip'j, k ,(x)} together constitute a biorthogonal 
wavelet basis o/L 2 (M). 
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• The discrete HT correspondences g'[k] = (d * g)[k] and g'[k] = (d * g)[k] 
hold. 

The above construction also exhibits the following properties: 

• The two biorthogonal systems have the same order, and the same Riesz 
bounds. 

• If the pair (cp, cp) satisfy the biorthogonality relation, then so do (cp', (p'). 
Indeed, using the identity TLf3" +l ^ 2 {x) = — A°^ 2 /3"(x), we can express the 
inner-product (cp', cp'(- — k)) as 

(0?+i/2 * Vo, (Pr+1/2 * <Po)(- - k)) 
= (W(^+l/2 * <Po),W(/Jf + i /2 * <Po)(- - k)} 
= (-A° /2 * cpo) , -A° /2 * cpo) (• - fc)) 
= (<P,<P(- - fc)), 
which establishes the assertion. 

• The lowpass filters on both the analysis and synthesis side are "half- 
sample" shifted versions of one another, and are related via the modulation 
of the discrete HT filter: 

H{z- X ) - Di-z-^H'iz- 1 ), 
H'(z) = D{-z- 1 )H{z). 

In particular, the filter is "half-sample" delayed on the analysis side, 
whereas on the synthesis side the filter has a "half-sample" advance. 

• The highpass filters on both the analysis and synthesis side are related 
through the FD filter as 

G{z- X ) = D(z)G'(z- 1 ) 1 
G'(z) = D(z)G(z). 

• If the analysis and synthesis filters of the original biorthogonal system 
satisfy the PR conditions 

Giz-^Giz) + H{z-^)H{z) = 1, 
G(z~ 1 )G(—z) + H{z- 1 )H{-z) = 0, 

then so do the filters of the HT pair. Indeed, since D(z)D(z~ 1 ) = 1, we 
have 

G'(z~ 1 )G'(z)+H'(z~ 1 )H'(z) 

= D(z- 1 )D(z)G(z- 1 )G'(z) 
+ D(-z)D(-z- 1 )H(z- 1 )H(z) 
= G(z- 1 )G(z) + Hiz'^Hiz). 
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Similarly, G' {z^G' {-z) + H' {z-^H' {-z) = 0. 

Note that above properties relate to a common theme: the unitary nature of 
the operators H and A involved in the wavelet and the filterbank construction, 
respectively. 

5 ID IMPLEMENTATION 

Signal Pre- filtering: In order to implement the DT-CWT, we need to employ 
two parallel wavelet decompositions corresponding to the wavelets ip(x) and 
ip'(x). Moreover, to have a coherent signal analysis — same input applied to 
both wavelet branches — we need to project the input signal f(x) separately onto 
V(<p) and V((p') before applying the respective DWTs. In particular, given a 
finite-energy input signal f(x), we consider its orthogonal projection fo{x) = 
co[k]<p(x — k) onto the space V(cp). The J-level wavelet decomposition of the 
signal fo{x) is subsequently given by 

]T djlktyjtix), (22) 

where the wavelet coefficients dj [k] , and the coarse approximation coefficients 
c.j[k] are recursively derived from the projection coefficients co[k] using Mallat's 
filterbank algorithm [22]. 

However, in practice one has access only to the discrete samples of the input 
signal f(x); let {/[fc]}feez be such (uniform) signal samples. It turns out that 
by assuming the input signal f{x) to bandlimitcd, a particularly simple digital 
filtering algorithm for computing the projection coefficients is obtained: 

co[k] = (f*p)[k], (23) 

where the frequency response P(e JaJ ) of the digital filter p[k] is uniquely specified 
by the restriction P(e- J ") = (p(w) for w G (— tt,tt) (derivation details in §11.3). 

As for the second branch, the input signal is projected onto the correspond- 
ing approximation space V(<p'): the same type of pre-filtering is applied with an 
appropriate modification of the frequency response, i.e., cp'(u>) is used instead 
of <p(w). To implement the filters for finite input signals, we use a FFT-based 
algorithm, similar to the one used in [3] for implementing the DWT filters. 

Analysis & Reconstruction: To simplify the notation, we shall henceforth 
use matrix notation to represent the linear transformations associated with the 
discrete DT-CWT. For instance, corresponding to an input signal f £ ~R N , the 
least-square projections are specified by Cq = Pf and c'q = P'f , where P and 
P' are the N x N circulant matrices corresponding to the two pre-filters. 

Let (h,g,h,g) and (h',g',h',g') be the set of perfect-reconstruction filters 
associated with the biorthogonal systems (cp,V>, (p,ip) and (cp', ip', cp', if/) re- 
spectively. The lowpass {(cf,c'f)} and the highpass {(cf* ,c'f)} subbands at 
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successive levels i — 1 , . . . , J are then given by the recursive interbank decom- 
positions: 

c'f = F & ,c'f_ 1 , c'f = F § ,c'f_i, (24) 

where and Fg (resp. F^, and F g >) denote the composition of the down- 
sampling matrix and the DWT matrix representing the lowpass and highpass 
analysis filters of the first (resp. second) channel. The complex wavelet sub- 
bands wi, . . . , wj are then specified by w t = cf + jc'f. In fact, the analysis 
can be summarized by the single frame operation 

T:f ~ {c L j,c' L J ,w 1 ,...,Wj). (25) 

from a lower-dimensional space to a higher-dimensional space: dim(Tf) > 
dim(f ). 

In several signal processing applications (e.g., denoising) one also needs to 
perform an inverse transform, that is, reconstruct the denoised signal from the 
processed complex wavelet coefficients. Since T is realized through the concate- 
nation of bases, it is injective: Tf = Tf only if / = /; however, as result of the 
redundancy, T exhibits non-unique left-inverses. In our case, we use a simple 
left-inverse: 

I* : (c L j, c' L j, wi, . . • , wj) ~ f = ^(P- 1 ^ + P'-V L ), 

where (cq ,c'q) are obtained via the recursion 

cf = F /l cf +1 + F 9 9te(w i+ i), 

c'f = F fc ,c'f +1 + F fl ,3m(w i+1 ), (26) 

for i = J — 1, . . . , (D\t(z) and 3m(z) denote the real and imaginary components 
of z respectively). Here F^ and F g (resp. Fh> and F g /) represent the composition 
of the DWT matrix corresponding to the lowpass and highpass synthesis filters of 
the first (resp. second) channel and the upsampling matrix. In short, the above 
inversion operation essentially amounts to inverting the two parallel transforms 
and averaging the inverses. 

Remark: The role played by the two projection filters P(e JW ) and P'^^) 
is critical as far as the issue of analyticity is concerned. Note that, while the 
analytic wavelet has an exact one-sided Fourier transform (by construction), 
the corresponding complex wavelet filter G(e^) +jG'(e : > UJ ) does not inherit this 
property naturally; it is only the combination of the projection and wavelet 
filters, Pa(e^) = P{e j ")G{e jw ) + jP'(&> u )G'(e> u ), that exhibits this property: 
Pa{& u ') — 0, for lo e (— 7r, 0]. Figure 1 shows the one-sided magnitude response 
of the filter. 
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Figure 1: Transfer function of the analytic wavelet filter P a (e : ' u ). 

6 GABOR-LIKE WAVELETS 

The "quantum law" for information — the principle that the joint time-frequency 
domain of signals is quantized, and that the joint time-frequency support of 
signals always exceed a certain minimal area — was enunciated in signal theory 
by Dennis Gabor [13]. He also identified the fact that the family of Gaussian- 
modulated complex exponentials (and their translates) provide the best trade-off 
in the sense of Heisenberg's uncertainty principle. 

The canonical Gabor transform analyzes a signal using the set of "optimally- 
localized" Gabor atoms: 

9m,nW = 7^ eX K-^ (X - TOT e ' nn( *~" lT) (2?) 

generated via the modulations-translations of a Gaussian-modulated complex 
exponential pulse [13, 1]. In particular, this paradigm involves the analysis of 
a finite-energy signal f(x) using the discrete sequence of projections c m ^ n = 
(fi9m,n) corresponding to different modulations and translations (m,n) E Z 2 . 

Note that the Gabor atoms have a fixed size (specified by the width T\ of 
the Gaussian window), and hence the associated transform essentially results 
in a "fixed-window" analysis of the signal. Moreover, the analysis functions 
{9m,n(x}} form a frame [21] and not a basis of L 2 (K); consequently, the recon- 
struction process involving the dual frame is often computationally expensive 
and/or unstable [1]. 
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6.1 Analytic Gabor-like Wavelets 



As a concrete application of the ideas developed in §4, we now construct a 
family of analytic spline wavelets that asymptotically converge to Gabor-like 
functions. In particular, we consider the family of semi-orthogonal B-spline 
wavelets that are better localized in space than their orthonormal counterparts, 
and that exhibit remarkable joint time-frequency localization properties [34] . 

In particular, consider the multiresolution in §2.2, generated by the fractional 
B-spline /3"(x). The transfer function of the wavelet filter that generates the 
so-called B-spline wavelet [35] associated with this multiresolution is specified 
by 

G^(e^) = e? u 'A a (-e> UJ )H?(-e- j0J ). (28) 

We denote the wavelet by ip"(x). The dual multiresolution (resp. wavelet) is 
specified by the unique dual-spline function $"(x) (resp. dual wavelet, denoted 
by C(x)). 

Following Corollary 4.3, it can be shown (proof provided in §11.2) that the 
family of B-spline wavelets ( x )}reR, of a fixed order a, and their duals are 
closed with respect to the HT: 

Proposition 6.1 (HT Pair of B-spline Wavelets) The HT of a B-spline (resp. 
dual-spline) wavelet is a B-spline (resp. dual-spline) wavelet of same order, but 
with a different shift: 

Hfc(x) = r T+l/2 {x). (29) 

The importance of this result is that it allows us to identify the analytic 
B-spline wavelet of degree a and shift r: 

*?(x)=W(x)+jW +1/2 (x). (30) 

In the sequel (cf. §7.3), we shall make particular use of this analytic spline 
wavelet. We would, however, like to highlight a different aspect: the remark- 
able fact that the wavelet *f>"(x) resembles the celebrated Gabor functions for 
sufficiently large a. Indeed, it was shown in [34] that the B-spline wavelets 
asymptotically converge to the real part of the Gabor function; by appropri- 
ately modifying the proof in [34], the following asymptotic convergence can 
established: 

^) ^ + oo Mcxp (- (a ^) cos - y — )■ < 31 > 

where M = 2M Q+1 Au /^27r(a + 1); a = V« + 1/Aw , with M = 0.670, lo q = 
—5.142 and Au> — 2.670. We recall that the asymptotic notation f a (x) ~ g a (x) 
signifies that f a {x)/g a {x) -^lasa^ +oo for all x. Immediately, we have the 
following result: 
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Proposition 6.2 (Gabor-like Wavelet) The complex B-spline wavelet *&"(x) 
resembles the Gabor function for sufficiently large a: 

*°(a;) ~ Mcxp (-^^r^) c j(-o*-^--t) _ (32) 

The above convergence happens quite rapidly. For instance, we have ob- 
served that the joint time-frequency resolution of the complex cubic B-spline 
wavelet (a = 3) is already within 3% of the limit specified by the uncertainty 
principle. Fig. 2 depicts the complex wavelets generated using HT pair of B- 
spline wavelets; the wavelets becomes more Gabor-like as the degree increases. 
Also shown in the figure is the magnitude envelope |^"(x)| of the complex 
wavelet which closely resembles the well-localized Gaussian window of the Ga- 
bor function. From a practical viewpoint, this means that one could use the 
non-redundant and numerically stable multiresolution spline transforms to ap- 
proximate the Gabor analysis. 

Remark: While the B-spline wavelets tend to be optimally localized in space, 
we have already observed that they are not orthogonal to their translates. The 
reconstruction therefore requires the use of some complementary dual functions. 
The flip side is that these dual-spline wavelets have a comparatively poor spatial 
localization, that deteriorates as the degree increases. This is evident in Fig. 3, 
which shows quadrature pairs of such wavelets of different degrees. However, 
we should emphasize that the dual (synthesis) wavelets have the same math- 
ematical rate of decay as their analysis counterpart, and that the associated 
reconstruction algorithm is fast and numerically stable. 



6.2 Gabor-like Transform 

The dual-tree Gabor-like transform is based on the analytic B-splinc wavelet 
$o(x) — ipo( x ) + jV'i/2( x )i where the degree a is sufficiently large (the choice 
r = is arbitrary). The analysis and synthesis filters for the first and second 
channel are as specified below [35]: 

• First channel: 

H(z) = 2- {a+1 \l + z)^{l + z- 1 )^, 
G(z) = zA a {-z)H{-z- 1 ), 
H(z) = H(z)A a (z)/A a (z 2 ), 

G(z) = G(z)/A a {z 2 )A a (-z). (33) 
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Real Part 

Imaginary Part 

Modulus 



(b) a = 6 



Figure 2: HT pairs of B-spline wavelets. In either case, Blue (solid line): iPq(x), 
Red (broken line): ip", 2 (x), Black (solid line): (x) + if;", 2 (x)\ 
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Figure 3: HT pairs of dual-spline wavelets. In either case, Blue (solid line): 
4>o(x), Red (broken line): tpy 2 (x). 

• Second channel: 

H'(z) = 2- { - a+1 \\ + z)t(l + ^- 1 )f+ 1 , 
G'(z) = zA a (~z)H'(-z~ 1 ), 
H'{z) = H'~(z)A a (z)/A a (z 2 ), 

G'(z) = G'{z)/A a {z 2 ) a A{-z). (34) 

The DT-CWT corresponding to this Gabor-like wavelet would then result 
in the analysis of the input signal f(x) in terms of the sequence of multiscale 
projections (/, v2"^&o (2 m • —It)) onto the (normalized) dilated-translated tem- 
plates of the Gabor-like wavelet (a;). Note that here the Gabor-like wavelet 
is used for analysis, whereas its dual is used for synthesis. The correspond- 
ing DWTs (F^,Fg,Fh,F g ) and (F^,,Fg/,Ffe/,F ff /) are efficiently implemented 
using a practical FFT-based algorithm, outlined in [3]. This method is exact 
despite the infinite support of the underlying wavelets, and achieves perfect- 
reconstruction up to a very high accuracy. The pre-filters, P and P', are also 
implemented in a similar fashion. An added advantage of the frequency domain 
implementation is that the execution time is independent of the order of the 
spatial filters. Moreover, the filters in (33) need to be pre-computed once and 
for all in order to apply the transform to different signals (of a fixed length) . 
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7 BIVARIATE EXTENSION 



Next, based on ideas similar to those of Kingsbury [27], we construct 2D com- 
plex wavelets, and 2D Gabor-likc wavelets in particular, using a tensor-product 
approach. Moreover, we also relate the real and imaginary components of the 
complex wavelets using a multi-dimensional extension of the HT. 

Separable Biorthogonal Wavelet Basis: Biorthogonal wavelet bases of L 2 (R) 
can be combined to construct a biorthogonal wavelet basis of L 2 (M 2 ). The 
underlying principle used to construct such a basis using tensor-products is as 
follows [21]: 

Theorem 7.1 Let (ip p ,^) p ) be the primal and dual wavelets of a biorthogonal 
wavelet basis o/L 2 (IR) ; with corresponding scaling functions (cp p ,cp p ). Simi- 
larly, let (ipq,ip q ) constitute another biorthogonal wavelet basis with correspond- 
ing scaling functions ((p q ,(p q ). Consider the following separable wavelets and 
their duals 

ip!(x) = (pp(x)ip q (y), ^i(x) = <p p {x)xi> q (y), 
ip 2 {x) = i/>p{x)<p q (y), tp 2 (x) = ^ p (x)(p q (y), 
tp 3 (x) = i/> p {x)ip q (y). tp 3 (x) = i[) p (x)%(y). (35) 

Then the dilation-translations of(ipi(x), ip2(x), ip3 (x)) and (^i(cc), tp2{x), ip 3 (x)) 
together constitute a biorthogonal wavelet basis o/L 2 (IR 2 ). 

The functions ipi(x), ip 2 ( x ) and ip 3 {x) are popularly referred to as the 'low-high' 
(LH), 'high- low' (HL) and 'high- high' (HH) wavelets, respectively, to emphasize 
the directions along which the lowpass scaling function and the highpass wavelet 
operate (here x = (x,y) denote the spatial coordinates). Note that the primal 
and dual approximation spaces for the above construction are V((p p ) (g> V((p q ) 
and V(<p p ) <8> V((p q ) respectively, where V((p p ) ® V((p q ) denotes the subspace 
span{cp p (- - m)q) q (- - n)} (m ,„ )eZ 2. 

7.1 Wavelet Construction 

A drawback of 2D separable wavelets is their preferential response to horizontal 
and vertical features. Fig. 4 shows the three separable wavelets arising from the 
separable construction. The pulsation of the LH and HL wavelets are oriented 
along the directions along which the constituent ID wavelets operate. However, 
the HH wavelet, with its constituent ID wavelets operating along orthogonal 
directions, does not exhibit orientation purely along one direction; instead it 
shows a checkerboard appearance with simultaneous pulsation along the diago- 
nal directions. 

This is exactly where the analytic wavelet tp a — ip + jH{ip}, with its one- 
sided frequency spectrum comes to the rescue: if instead of employing separable 
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wavelets of the form i))(x)ip(y), complex wavelets of the form ip a (x)tp a (y) are 
used, then the corresponding spectrum -ip a (uj x )ip a (ujy) will have only one pass- 
band, and consequently the real wavelets 9iz('ip a (x)'ip a (y)) and 3m(ip a (x)ip a (y)) 
will indeed be oriented. 

The motivation then is to use HT pairs of ID biorthogonal wavelets to con- 
struct oriented 2D wavelets. In particular, we do so by appropriately combin- 
ing four separable biorthogonal wavelet bases using Theorem (7.1). To begin 
with, we immediately identify the two scaling functions (p p (x) = (p(x) and 
(p q (x) = cp'(x), associated with the analytic wavelet ip a {x) = "4>(x) + jip'(x), 
where i/>'(x) ~ Tiip(x). This naturally leads to the possibility of four separable 
biorthogonal wavelet bases corresponding to the following possible choices of 
approximation spaces: V((p) ® V(q>), V(q>) <£> l^(cp'), V((p') £g> V((p) and V((p')<£> 
V((p'). In fact, as will be demonstrated shortly, we will employ all of these to 
obtain a balanced construction. 

First, we identify the separable wavelets corresponding to the four scaling 
spaces: 

i>i{x) = ip{x)ip{y), ip 4 {x) = (p(x)ip'(y), 
ip 2 {x) = ip(x)(p(y), ip 5 (x) = i)(x)<p'(y), 
ip 3 (x) = ip(x)i>(y), ip 6 (x) = i)(x)ip'(y), 



ip 7 (x) = (p'(x)tp(y), i>w(x) = (p'(x)ip'(y), 
tp$(x) = ip'(x)<p(y), i>u{x) = i>'(x)<p'(y), 
Mx) = 1/(x)il>(y), Vi2(x) = i//(xW(y). (36) 

The corresponding dual wavelets ^i(x) are specified identically except that the 
dual wavelets are used instead of the primal ones. Finally, by judiciously using 
the one-sided spectrum of the analytic wavelet ip a (x) = tp(x) + jip'(x), and 
by combining the four separable wavelet bases (36), we arrive at the following 
wavelet specifications: 

V 2 (x) = ip a {x)<p'{y) = ip 5 {x) +ji/> 11 (x), 
* 3 (je) = (p{x)ip a (y) = i>i(x)+ji/)4(x), 

^ 4 (x) = (p'(x)tpa(y) = ll>7{x)+ji> 10 {x), 

v 2 (37) 
= l V2 ) +J { V2 J' 

( ip3(x) + ipi 2 (x) \ . ( tpejx) - ^ 9 (x) \ 
V V2 J +J \ V2 J' 
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Figure 4: Wavelets associated with the separable basis. The figure shows the 
LH, HL and HH wavelets in the space domain. 

The dual complex wavelets, ^(a;), are specified in an identical fashion using 
the dual wavelets iptix). Importantly, the above construction is complete in 
the sense that it involves all the 4x3= 12 separable wavelets of the four 
parallel multircsolutions. The factor 1 /v2 ensures normalization: the real and 
imaginary components of the six complex wavelets have the same norm. 

7.2 Directional Selectivity and Shift-Invariance 

A real wavelet has a bandpass spectrum that is symmetric w.r.t. to the origin. 
As a result, for the wavelet to be oriented, it is necessary that its spectrum 
be bandpass only along one preferential direction. We claim that the real and 
imaginary components of the above complex wavelets are oriented along the 
primal directions 9\ = 9 2 = 0, 9 3 = 9 4 = ir/2,9 5 = tt/4, and 9 6 — 37r/4, 
respectively. Indeed, it is easily seen that the support of ^i(uj) and Vf^Cw) is 
restricted to the half-plane {(u) Xl LJy) : lj x > 0}, since their Fourier transform 

can be written as ^(w) = (1 + sign(w a; )) £He(^ , fe)(a;). As it is necessary for 
the real functions 9 = te(\E r fc) and 3m(*i>k) to have symmetric passbands, the claim 
about their orientation along the horizontal direction then follows immediately. 
The orientation of the components of the wavelets 'J 3 (x) and ^4(x) along the 
vertical direction follows from a similar argument. 

As far as the wavelet ^(te) is concerned, note that ^s(uj) = (1 + signf^)) (1 + sign(wj,)) ip(ui x )'il)(}jJy) ■ 
As a consequence, the support of Vf^oj) is restricted to the quadrant {(lo Xi uj y ) : 
uj x > 0,ujy > 0}. The symmetry requirements on the spectrums of V\c(^z)(x) 
and Dtm(^ , 5)(a;) then establish their orientation along ir/4. A similar argument 
establishes the orientation of the real components D\c(^a)(x) and 3m(^Q)(x) 
along 37r/4. 

The above-mentioned directional properties allude to some kind of analytic 
characterization of the complex wavelets. Indeed, akin to the ID counterpart, 
it turns out that the components of the above complex wavelets can also be re- 
lated via a multi-dimensional extension of the HT that provides further insights 
into the orientations of the wavelets. In particular, we consider the following 
directional version of the HT [14]: 

H e f(x) ^ -jsign(o; T Me )/(a;), (38) 
specified by the unit vector ug = (cos 9, sin 9) pointing in the direction ^ 
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Figure 5: 2D Gabor-like Wavelets. Left: Real component of the six complex 
wavelets, Right: Magnitude envelope of the six complex wavelets. The diago- 
nally placed wavelets are identical, they are used twice to balance the represen- 
tation. 

9 < 7r. That is, the directional HT is performed with respect to the half-spaces 
{uj : u> T ug > 0} and {u> : u> T ug < 0} specified by the vector ug, and it maps 
the directional cosine cos(u>^a;) into the directional sine sm(ujx). Based on the 
wavelet definitions (37), the following correspondences (for a proof see §11.4) 
can then be derived: 

Proposition 7.2 The real and imaginary components of the complex wavelets 
^k( x ) form directional HT pairs. In particular: 



A significant problem with the decimated DWT is that the critical down- 
sampling makes it shift- variant. The redundancy of the dual-tree transform has 
been successfully exploited for partially mitigating this shift-variance problem 
[17, 19]. Our design further mitigates this shift- variance problem by using a 
finer sub-sampling scheme in the 0° and 90° directions. Observe that ^i(x) as 
* 2 (x, y ~ 1/2) owing to the fact that 0?{x) w f)" J2 (x + 1/2). These wavelets 
provide a finer sampling in the y-direction. Similarly, the vertical wavelets give 
us a finer sampling in the x-direction. 

7.3 Gabor-like Wavelets 

Daugman generalized the Gabor function to the following 2D form 



Jm(tf fc (x)) = Hg k 9\t(^k(x)), 1 < k ^ 6. 



(39) 



Sf(x) 
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(40) 
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involving the modulation of an elliptic Gaussian using a directional plane-wave, 
to model the receptive fields of the orientation-selective simple cells in the visual 
cortex [9]. 

We are particularly interested in the dual-tree wavelets, denoted by &i(x; a, r), . . . , ^ ( x ; a > T )> 
derived from the quadrature B- spline wavelets tl> p (x) = i>"(x) and tp q (x) = 
^t+i/2( x )- These complex wavelets inherit the asymptotic properties of the 
constituent spline functions. Indeed, by appropriately modifying the proof in 
[34], it can be shown that 

ff(l) ~7^" p (-^) < 41 > 

for sufficiently large a, where a = \Ja + I/2\/3. This, combined with (32), then 
results in the following asymptotic characterization: 

Proposition 7.3 f2D Gabor-like Waveletsj The complex wavelets ^k(x;a,r) 
resemble the 2D Gabor functions for sufficiently large a: 



&i(x; a,r) — Mie 



@2{x; a, t) ~ Mie 



^3 (x; a, t) <~ Mie 



&4(x; a, t) <~ Mie 



^(x; a, t) ~ M 2 e 



%?()(x; a, t) <~ M 2 e 

where M x = 2\/3M a+1 A a>o/7r(a + 1); M 2 = 2M 2( " +1) Awg/7r(a + 1); <j x = 
yja+ l/Awo; and cr 2 = \/ (a + l)/6. 

We call the wavelets "Gabor-like" since they form approximates of 2D Gabor 
functions similar to the ones proposed by Daugman (40). The dual-tree trans- 
form (cf. §8) corresponding to a specific family of such Gabor-like wavelets 
(fixed a and r) results in a multiresolution, directional analysis of the input im- 
age f(x) in terms of the sequence of projections (/, 2 %c S^(2 % x — m; a, r)). Fig. 5 
shows the 2D Gabor wavelets corresponding to a = 6 and r — 0. The ensemble 
shows the modulus |S4(x;6,0)| and the real component *He(%(a;; 6, 0)) of the 
six complex wavelets; the former shows the pulsations of the directional plane 
waves, whereas the latter shows the elliptical Gaussian envelopes. 



/ (s-l/2) 2 (y-T) 2 \ 

('-1/2) J (»-r-l/2) 2 \ 

V "I °2 J gj(w x--£—nT) 
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7.4 Discussion 



Before moving on to the implementation, we digress briefly to discuss certain 
key aspects of our construction: 

• Directionality: The six complex wavelets in Kingsbury's DT-CWT scheme 
are oriented along the directions: ±15°, ±45°, and ± 75° [19]. Though we 
use similar separable building blocks in our approach, our wavelets are 
oriented along the four principal directions: 0,7r/4,7r/2 and 3ir/4. The 
added redundancy along the horizontal and vertical directions yields better 
shift-invariance along these directions. Alternatively, we could also have 
applied Kinsbury's construction to obtain Gabor-like wavelets orientated 
along ±15°, ±45°, and ±75°. 

• Localization Vs. Frame Bounds: In this paper, we placed emphasis on 
time-frequency localization, and were able to construct new wavelets that 
converge to Gabor-like functions. These basis functions should prove use- 
ful for image analysis tasks such as extraction of AM-FM information 
and texture analysis. However, the price to pay for this improved local- 
ization is that the associated transform — in contrast with the transforms 
constructed by Kingsbury et al. [27, 18] is no longer tight, and con- 
sequently requires a different set of reconstruction filters. Nevertheless, 
the tightness of the frame-bounds — a desirable property for image pro- 
cessing applications such as denoising and compression — can, in principle, 
also be achieved within our proposed framework by replacing the B-spline 
wavelets with the orthonormal ones (Battle-Lemarie wavelets). 

• Analytic Properties: Our method of construction takes a primary wavelet 
transform and obtains an exact HT pair using a simple unitary mapping. 
The consequence is that all fundamental approximation-theoretic proper- 
ties of continuous-domain wavelets, such as vanishing moments and regu- 
larity, are automatically preserved, and that the associated filters inherit 
an exact one-sided response. We also obtain an explicit space-domain 
expression for the Gabor-like wavelets. 

• Multidimensional HT properties: The directional HT correspondences 
(39) for our complex wavelets follows as a direct consequence of the tensor- 
product construction. We would however like to note that there exist other 
multidimensional extensions of the HT as well: the "single-orthant" exten- 
sion of Hahn [15] involving the boundary distribution of analytic functions; 
the "hypercomplex" extension due to Billow et al. [6]; the "monogenic" 
signal due to Fclsberg et al. [11]; and the spiral- phase quadrature trans- 
form of Larkin et al. [20] . The last two in the list are closely related to the 
Riesz transform of classical harmonic analysis [32] . Design of directional 
wavelets based on these and other alternative extensions are a promising 
topic of research [7, 24]. 
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Figure 6: Block Diagram of the 2D Complex Wavelet Transform. 

8 2D IMPLEMENTATION 

Pre- filtering: The input signal has to be projected onto each of the four sep- 
arable spaces of the form V(q> p ) <g> V((p q ) before initiating the multiresolution 
decompositions. As in the ID setting, the orthogonal projection is achieved in a 
separable fashion using an appropriate prc-filtcr along each dimension. In par- 
ticular, if {f[k]}kez 2 be the uniform samples of a bandlimited input signal f(x), 
then the projection coefficients are given by c$ L [k] = (/ *p)[fe], where the sep- 
arable pre-filter p[k] is specified by 5Zp[fci, kJ\e~^ klU}x+k2UJ ^ = q>*(u x )<Ps{<* } y) 
for u> = (oj x ,u! y ) in (— n,ir) 2 . 

In general, there would be four such projections Cy L (n) = f *p n , 1 n ^ 4, 
corresponding to the 2D pre-filters p\,...,p± associated with the four approx- 
imation spaces. Note that the filters can be implemented efficiently through 
successive ID filtering along cither dimension. 

Analysis: We consider the implementation aspects for a finite input signal 
f e R MxN . The transform, corresponding to the complex wavelets (37), involves 
four separable DWTs with different filters applied along the x and y directions 
(cf. Table 8 for the list of filters), and result in four subbands at each decom- 
position level. Specifically, let cf L (n), c\ H {n), c^ L (n) and cf H (n),l ^ n 4, 
denote the low-low, low-high, high-low and high-high subbands, respectively, of 
the four DWT decompositions at resolution i = 1 , . . . , J. The low-low subbands 
Cq L (u) arc identified as the four set of pre-filtered signals cf L (n) = P„f, with 
P„ being the (block) circulant matrices associated with the 2D pre-filters. The 
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coarser subbands at levels 



. , J are then given by 



cf L (n) = F n (h x ,h y )cf^(n) 

cf H (n) = F n (h x ,g y )cf^(n) 

cf L (n) = F n {g x ,h y )Ci-M 

cf H (n) = F n (fe,ff B )cfii(n), (43) 

where F n (q x ,q y ) denotes the composition of the nth DWT matrix (employing 
analysis filters q x and q y in the ^-direction and y-direction) , and the downsam- 
pling matrix. 

The complex subbands w» = (wj, . . . , wf), 1 < i < J, are specified by w, = 
A«rC» + JA3&, where 

Ct = (c^(l),c^(2),c^(l),c^(3),cr ( l),cf«(4)), 
6 = (c^(3),c^(4),c^(2),c^(4),cr ( 2),cf«(3)), 

are obtained through a particular permutation of the 12 highpass subbands; and 
the block matrices and A3 are specified as 
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(45) 



In short, the transform can be formally summarized via the frame operation 

T:f ~(c J (l),...,cj(4),wi,...,w J ) (46) 

involving the sequence of transformations: projections Pi, ... , P4; discrete wavelet 
transforms F l7 ...,F4; permutation II, and orthonormal transformations A^ 
and A3. Figure 6 provides a schematic of these sequence of transformations. 



Reconstruction: Note that the permutation 

II : {cf*(n),c?>),c?*(n)} Kn< 4 - (Ci.Ci) (47) 

involved in (44) is invcrtiblc, and that the matrices A<h and A3 are orthonormal, 
with corresponding inverses given by A^ and A3 respectively. Starting with the 
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Figure 7: Directional decomposition (one-level) of a synthetic image (Octagon) 
and a natural image (Cameraman) using the Gabor-like transform. Ordering 
of the subbands in either case: First column : |w 1 |,|w 2 |; Second column : 
|w 3 |, |w4|; Third column : |w 5 | and |w 6 |. 
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complex wavelet subbands Wi, . . . , wj, the highpass subbands c*(l), . . . , c^(4) 
corresponding to the bands b — HL, LH, and HH, are then computed from 
the vectors Ci = A^£Re(wi) and & — Aij3m(wj), via the permutation II -1 at 
levels i — 1, . . . , J. These, along with the lowpass subbands Cj L (l), . . . , Cj L (4), 
are then used to reconstruct the projected signals Cq L (1), . . . , Cq L (4) using the 
recursion 

c^ L (n) = F„(^,^)cf + L 1 (n) +F n (h x ,g y )c^ 1 (n), 

+F n (g x , h y )cfji („) + F n (g x , g y )c%[ (n) (48) 

for i = J — 1, ... ,0. Here, F n (m x ,m y ) represents the composition of the up- 
sampling matrix and the synthesis matrix corresponding to the nth DWT, 
with filters m x and m y in the x-direction and y-direction, respectively, as 
specified in Table 8. The input signal samples are finally recovered as f = 

2D Gabor-like Transform: The Gabor-like transform is based on the ana- 
lytic B-splinc wavelets specified in §7.3, where the complex subbands wf [ra] 
represent the directional decompositions of the input image along the four pri- 
mal directions using the optimally-localized Gabor-like wavelets a, r) at 
different resolutions. The filtcrbank analysis (43) and synthesis (48) operations 
arc implemented in a separable fashion using the ID spline DWT filters specified 
in (33) and (34). 

Fig. 7 shows the magnitude response of the six complex wavelet subands 
obtained by applying our Gabor-like transform to a synthetic and a natural im- 
age. In particular, the wavelet subbands corresponding to the synthetic image, 
with directional edges along 0,7r/4,7r/2 and 37r/4, highlight the directional- 
selectivity of the transform. The simulation was carried out in MATLAB 7.5 on 
a Macintosh 2.66 GHz Intel dual-core system. The average execution time for 
one-level wavelet analysis and reconstruction (including pre- and post-filtering) 
of a512 x 512 image is 1.2 seconds, and the reconstruction error is of the order 
of lO" 16 . 
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Figure 8: Analysis and synthesis filters corresponding to the four multiresolu- 
tions. 
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9 CONCLUDING REMARKS 



The primary objective of this contribution was to combine the attractive fea- 
tures of Gabor analyses and multiresolution wavelet transforms into a single 
theoretical framework, and to provide a fast algorithm for the same. Specifi- 
cally, we proposed a formalism for constructing exact HT pairs of biorthogonal 
wavelets based on (i) the B-splinc factorization theorem, and (ii) a natural dis- 
cretization of the continuous HT filter identified via the action of the HT on 
fractional B-splines. Based on this methodology, analytic wavelets resembling 
the Gabor function were then designed using HT pair of B-spline wavelets. 

We then extended our scheme to 2D: starting from HT pair of ID biorthogo- 
nal wavelet basis, we constructed directional complex wavelets by appropriately 
combining four separable biorthogonal wavelet bases. In particular, we related 
the real and imaginary components of the complex wavelets using a directional 
extension of the HT. The particular family of wavelets constructed using B- 
splines was shown to resemble the directional Gabor function family proposed 
by Daugman. Finally, we demonstrated how the discrete Gabor-like transforms 
could be implemented using fast FFT-based filtcrbank algorithms. 
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11 APPENDIX 

11.1 Proof of Theorem 4.2: 

We begin with the following sequence of equivalences 

Hip(x/2) = g[k]H(p(x - k) 

feez 

= £ff[fc](W*<Po)(a-*) 

feez 

= EsM ( E d N(/?" + i/2 * <Po)0 -n)\x-k) 
feez Vnez / 

= ^2(9*d)[m](p'(x-m), 

mGZ 

based on (12), and the linearity, associativity, and commutativity of the under- 
lying convolution operators. The sufficiency part of the theorem then follows 
immediately: if g'[k] = (g * d)[k], then tp'(x) = Htp(x). 
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Conversely, let tp'(x) — Hip(x), so that ^g'[k]<Q'(x — k) = ^2(g*d)[k](p'(x — 
k). Now, since {cp'(- — n)} forms a Riesz basis of the subspace V((p') = 
span^ 2 {cp'(- — k)}k£Z, every element in V(q>') necessarily has a unique rep- 
resentation. Hence, g' [k] — (g * d) [k] . 

11.2 Proof of Proposition 6.1: 

The primal scaling functions can be trivially factorized: <p(x) = (/3"*cp )(x) and 
cp'(x) = (/3" +1 / 2 * <Po)(x), where cp is the Dirac delta distribution. Similarly, 
the dual scaling functions can be factorized as <p(x) = * cpo)(a:) and cp'(x) = 



(/ 3 r+i/2 * <Po)(x), where cp = £g a [fc]<5(- ~ k ) with E q a [k]e~ jwk = l/A a (e juJ ). 



Note that in the latter case we have particularly used the fact that A a 
and hence q a [k], are independent of r. 

The proposition then follows from Corollary (4.3) since the wavelet filters 
satisfy the sufficiency conditions: g'[k] = (d* g)[k] and g'[k] = (d* g)[k], respec- 
tively. Indeed, from (14) and (28), we have 



The other condition G'{c^) = D{e 3u> )G{e ju> ) can be similarly derived. 
11.3 Derivation of Equation (23): 



gives the orthogonal projection of f(x) onto V(q>). The solution to the above 
problem is explicitly given by Pv( v )f(x) = ^2c a [k](p(x — k) where the coeffi- 
cients are specified by Co[k] = (/, cp(- — k)). Here cp(x) denotes the dual of <p(x) 
that satisfies the biorthogonality criterion (cp, (p(- — n)) = S[n). Moreover, under 
the constraint that cp(a;) e V((p), we recover a unique dual that is specified by 
the Fourier transform (p(w) = $(w)/E l$( w + 27r/c)| 2 [33]. 

Next, using the Poisson summation formula, we derive the expression Co(e 3UJ ) = 

Ekgz (/ * ( P T )( W + 27rn) for the (discrete) Fourier transform of Co[k]. The ban- 
dlimited model f(x) — £ fW\ sinc(a; — k) finally results in the simplification 



G'{d w ) = e ju 'A a (-e : >")H? +1/2 (-e- j ") 

= c^A a (-c^)D(e juJ )H?(-c- juJ ) 
= L>(c^)G(e^). 



It is well-known that the least-square approximation operator Pv(<p) ■ L 2 (M) 
V((p), defined by 

p v(<p)f = ar g . min 11/ - /oil (< 



(49) 



Co(e^) 




= F(e^)P(e^) 



(50) 
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where P(e J ") equals (p(w) on and F(e^ UJ ) is the Fourier transform of 

/[*]■ 



11.4 Proof of Proposition 7.2: 

We establish the correspondence for the wavelets ^i(x) and ^(a;) (the rest can 
be derived similarly) . The correspondence for the former is direct: Ho'$lc('i'i(x)) = 
HoMx)}<p(y) = ip'{x)<p{y) = Jm(*i(aj)). 

Next, note that the Fourier transforms of fRe(^f^(x)) and Um^^x)) can be 
written as 

9te(* 5 )(w) = + sign(w 2 )sign(c^))V>(w K )^(u;, y ), and 

V2 

3m(* 5 )(o>) = --^=(sign(o; x ) + sign(w y ))^(w x )V;(wj / ). (51) 

The correspondence 3m( 1 $5(x)) = H^/^e^^^x)) then follows from the iden- 
tity (sign(w x ) + sign(^)) = sign(w a: + w„)(l + sign(w a; )sign(^)) . 



References 

[1] M. J. Bastiaans, Gabor's expansion of a signal into Gaussian elementary 
signals, Proc. IEEE 68 (1980), 538-539. 

[2] J. J. Benedetto, Harmonic Analysis and Applications, CRC Press, 1996. 

[3] T. Blu and M. Unser, The fractional spline wavelet transform: Definition 
and implementation, Proc. IEEE International Conference on Acoustics, 
Speech, and Signal Processing (ICASSP'00) (2000), 512-515. 

[4] , A complete family of scaling functions: The (a, r) -fractional 

splines, Proc. IEEE International Conference on Acoustics, Speech, and 
Signal Processing (ICASSP'03) VI (2003), 421-424. 

[5] R. Bracewcll, The Fourier Transform and Its Applications, McGraw-Hill, 
1986. 

[6] T. Biilow and G. Sommcr, Hypercomplex signals - A novel extension of the 
analytic signal to the multidimensional case, IEEE Trans. Signal Process. 
49(11) (2001), 2844-2852. 

[7] W. L. Chan, H. Choi, and R.G. Baraniuk, Coherent multiscale image pro- 
cessing using dual-tree quaternion wavelets, IEEE Trans. Image Process. 
17 (2008), 1069-1082. 

[8] C. Chaux, L. Duval, and J.C. Pcsquet, Image analysis using a dual-tree 
M-band wavelet transform, IEEE Trans. Image Process. 15 (2006), no. 8, 
2397-2412. 



33 



[9] J. G. Daugman, Two-dimensional spectral analysis of cortical receptive field 
profile, Vision Research 20 (1980), 847-856. 

[10] P. F. C. de Rivaz and N. G. Kingsbury, Bayesian image deconvolution and 
denoising using complex wavelets, Proc. IEEE International Conference on 
Image Processing (ICIP'01) 2 (2001), 273-276. 

[11] M. Felsberg and G. Sommer, The monogenic signals, IEEE Trans. Signal 
Process. 49(12) (2001), 3136-3144. 

[12] F. C. A. Fernandes, R. L. C. Van Spaendonck, and C. S. Burrus, A new 
framework for complex wavelet transforms, IEEE Trans. Signal Process. 51 
(2003), no. 7, 1825-1837, 43. 

[13] D. Gabor, Theory of communication, J. Inst. Elect. Eng. 93 (1946), 429- 
457. 

[14] G. H. Granlund and H. Knutsson, Signal Processing for Computer Vision, 
ch. 4, Dordrecht, The Netherlands: Kluwer, 1995. 

[15] S.L. Hahn, Multidimensional complex signals with single- orthant spectra, 
Proc. IEEE 80(8) (1992), 1287-1300. 

[16] S. Hatipoglu, S.K. Mitra, and N.G. Kingsbury, Image texture description 
using complex wavelet transform, Proc. IEEE International Conference on 
Image Processing (ICIP'00) 2 (2000), 530-533. 

[17] N. G. Kingsbury, Shift invariant properties of the dual-tree complex wavelet 
transform, Proc. IEEE International Conference on Acoustics, Speech, and 
Signal Processing (ICASSP'99) (1999). 

[18] , A dual-tree complex wavelet transform with improved orthogonality 

and symmetry properties, Proc. IEEE International Conference on Image 
Processing (ICIP'06) 2 (2000), 375-378. 

[19] , Complex wavelets for shift invariant analysis and filtering of sig- 
nals, Journal of Applied and Computational Harmonic Analysis 10 (2001), 
no. 3, 234-253. 

[20] K. G. Larkin, D. J. Bone, and M. A. Oldficld, Natural demodulation of 
two-dimensional fringe patterns. I. General background of the spiral phase 
quadrature transforms, J. Opt. Soc. Am. A 18(8) (2001), 1862-1870. 

[21] S. Mallat, A Wavelet Tour of Signal Processing, San Diego, CA: Academic 
Press, 1998. 

[22] S. G. Mallat, A theory for multiresolution signal decomposition: The 
wavelet representation, IEEE Trans. Pattern Anal. Mach. Intcll. 11 (1989), 
no. 7, 674-693. 

[23] Y. Meyer, Ondelettes et operateurs i: Ondelettes, Hermann, Paris, 1990. 



34 



[24] S.C. Olhedc and G. Metikas, The hyperanalytic wavelet transform, Imperial 
College Statistics Section Technical Report TR-06-02 (2008), 1-49. 

[25] M. S. Pattichis and A. C. Bovik, Analyzing image structure by multidimen- 
sional frequency modulation, IEEE Trans. Pattern Anal. Mach. Intcll. 29 
(2007), no. 5, 753-766. 

[26] I. W. Selesnick, The design of approximate hilbert transform pairs of wavelet 
bases, IEEE Trans. Signal Process. 50 (2002), no. 2, 1143-1152. 

[27] I. W. Selesnick, R. G. Baraniuk, and N. C. Kingsbury, The dual-tree com- 
plex wavelet transform, IEEE Signal Process. Mag. 22 (2005), no. 6, 123- 
151. 

[28] I.W. Selesnick, Hilbert transform pairs of wavelet bases, IEEE Signal Pro- 
cess. Lett. 8 (2001), no. 6, 170-173. 

[29] L. Sendur and I. W. Selesnick, Bivariate shrinkage with local variance es- 
timation, IEEE Signal Process. Lett. 9 (2002), no. 12, 438-441. 

[30] E. M. Stein, Singular Integrals and Differentiability Property of Functions, 
Princeton University Press, 1970. 

[31] E. M. Stein and R. Shakarchi, Complex analysis, Princeton University 
Press, 2003. 

[32] E. M. Stein and G. Weiss, Fourier Analysis on Euclidean Spaces, Princeton 
University Press, 1971. 

[33] M. Unser, Sampling— 50 years after Shannon, Proc. IEEE 88 (2000), no. 4, 
569-587. 

[34] M. Unser, A. Aldroubi, and M. Eden, On the asymptotic convergence of 
B-spline wavelets to Gabor functions, IEEE Trans. Inf. Theory 38 (1992), 
no. 2, 864-872. 

[35] M. Unser and T. Blu, Construction of fractional spline wavelet bases, Proc. 
SPIE Conference on Mathematical Imaging: Wavelet Applications in Sig- 
nal and Image Processing VII 3813 (1999), 422-431. 

[36] M. Unser and T. Blu, Fractional splines and wavelets, SIAM Review 42 
(2000), no. 1, 43-67. 

[37] , Wavelet theory demystified, IEEE Trans. Signal Process. 51 (2003), 

no. 2, 470-483. 

[38] B. Van der Pol, The fundamental principles of frequency modulation, Jour- 
nal IEE 93 (1946), 153-158. 

[39] J. Ville, Theorie et application de la notion de signal analytique, Cables 
and Transmissions 93 (1948), no. Ill, 153-158. 



35 



[40] R. Yu and H. Ozkaramanli, Hilbert transform pairs of biorthogonal wavelet 
bases, IEEE Trans. Signal Process. 54 (2006), 2119-2125. 



36 



